home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsiii / accform.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-12-13  |  7.3 KB  |  260 lines
  1. /*
  2. ACCURATE FORM-FACTOR COMPUTATION
  3.  
  4. Filippo Tampieri
  5. Cornell University
  6. */
  7. /*
  8.     The routines in this module provide a way of computing the
  9.     radiosity contribution of a source patch to a point on a
  10.     receiving surface.  The source must be a convex quadrilateral
  11.     and have a uniform radiosity distribution (or approximately
  12.     uniform).  For simplicity and clarity, the radiosity is
  13.     computed for a single wavelength.
  14. */
  15.  
  16. #include <stdio.h>
  17. #include <math.h>
  18.  
  19. static float computeFormFactor() ;
  20. static float computeUnoccludedFormFactor() ;
  21. static void splitQuad() ;
  22. static float quadArea() ;
  23.  
  24. #define X    0
  25. #define Y    1
  26. #define Z    2
  27.  
  28. /* set vector v to zero */
  29. #define VZERO(v)        (v[X] = v[Y] = v[Z] = 0.0)
  30. /* increment vector v by vector dv */
  31. #define VINCR(v, dv)    (v[X] += dv[X], v[Y] += dv[Y], v[Z] += dv[Z])
  32. /* scale vector v by a constant c */
  33. #define VSCALE(v, c)    (v[X] *= c, v[Y] *= c, v[Z] *= c)
  34. /* copy vector u to vector v */
  35. #define VCOPY(v, u)     (v[X] = u[X], v[Y] = u[Y], v[Z] = u[Z])
  36. /* compute the dot product of vectors u and v */
  37. #define VDOT(u, v)      (u[X] * v[X] + u[Y] * v[Y] + u[Z] * v[Z])
  38. /* compute the cross product of vectors s and t */
  39. #define VCROSS(v, s, t) (v[X] = s[Y] * t[Z] - s[Z] * t[Y], \
  40.                          v[Y] = s[Z] * t[X] - s[X] * t[Z], \
  41.                          v[Z] = s[X] * t[Y] - s[Y] * t[X])
  42. /* compute the length of vector v */
  43. #define VNORM(v)        (sqrt(VDOT(v, v)))
  44. /* subtract vector t from vector s and assign the result to v */
  45. #define VSUB(v, s, t)   (v[X] = s[X] - t[X], \
  46.                          v[Y] = s[Y] - t[Y], \
  47.                          v[Z] = s[Z] - t[Z])
  48. /* compute the average of vectors s and t */
  49. #define VAVG2(v, s, t)  (v[X] = 0.5 * (s[X] + t[X]), \
  50.                          v[Y] = 0.5 * (s[Y] + t[Y]), \
  51.                          v[Z] = 0.5 * (s[Z] + t[Z]))
  52.  
  53. typedef float Vector[3];
  54. typedef Vector Point;
  55. typedef Vector Normal;
  56.  
  57. /* represent a quadrilateral as an array of four points */
  58. typedef Point Quad[4];
  59.  
  60. /*
  61.     computeVisibility(Pr, S)
  62.     Point Pr;
  63.     Quad S;
  64.  
  65.     computeVisibility--takes as input a receiving point 'Pr' and a
  66.     source 'S' and estimates the visibility of 'S' as seen from 'Pr'.
  67.     It returns 0 for total occlusion, 1 for total visibility, and a
  68.     number representing the fraction of 'S' that can be seen from
  69.     'Pr' in the case of partial occlusion.
  70. */
  71. extern float computeVisibility();
  72.  
  73. static float eps_F, eps_A;
  74.  
  75. /*
  76.     computeContribution--takes as input the descriptions of a source
  77.     patch and a receiving point and returns the radiosity
  78.     contribution of the source to the receiver. The source is
  79.     described as a quadrilateral 'S', its UNIT normal 'Ns', and its
  80.     average unshot radiosity 'Bs'. The receiving point is described
  81.     by its location 'Pr', its UNIT normal 'Nr', and its reflectivity
  82.     'rho'. The parameter 'eps_B' controls the accuracy of the
  83.     computation in the presence of partial occlusion; smaller values
  84.     of 'eps_B' will result in more accurate estimates at the cost of
  85.     computation speed. The parameter 'minA' provides the means to
  86.     control the recursive subdivision of the source: no region of
  87.     area less than or equal to 'minA' will be subdivided.
  88. */
  89. float computeContribution(Pr, Nr, rho, S, Ns, Bs, eps_B, minA)
  90. Point Pr;
  91. Normal Nr, Ns;
  92. float rho, Bs, eps_B, minA;
  93. Quad S;
  94. {
  95.     Vector v;
  96.     float computeFormFactor();
  97.  
  98.     if(VDOT(Nr, Ns) >= 0.0)
  99.         /* the receiving point is oriented away from the source */
  100.         return 0.0;
  101.     
  102.     VSUB(v, Pr, S[0]);
  103.     if(VDOT(Ns, v) <= 0.0)
  104.         /* the receiving point is behind the source */
  105.         return 0.0;
  106.  
  107.     eps_F = eps_B / (rho * Bs);
  108.  
  109.     eps_A = minA;
  110.  
  111.     return(rho * Bs * computeFormFactor(Pr, Nr, S, Ns));
  112. }
  113.  
  114. /*
  115.     computeFormFactor--takes as input the descriptions of a source
  116.     patch and a receiving point and returns the form-factor from a
  117.     differential area centered around the receiving point to the
  118.     finite area source.  The source is described as a quadrilateral
  119.     'S' and its UNIT normal 'Ns'.  The receiving point is described
  120.     by its location 'Pr' and its UNIT normal 'Nr'. This routine
  121.     relies on the external function 'computeVisibility' to determine
  122.     a visibility factor.
  123. */
  124. static float computeFormFactor(Pr, Nr, S, Ns)
  125. Point Pr;
  126. Normal Nr, Ns;
  127. Quad S;
  128. {
  129.     int j;
  130.     Quad dS[4];
  131.     void splitQuad();
  132.     float Frs, vis, computeUnoccludedFormFactor(), quadArea();
  133.  
  134.     vis = computeVisibility(Pr, S);
  135.     if(vis <= 0.0)
  136.         Frs = 0.0;
  137.     else if(vis >= 1.0)
  138.         Frs = computeUnoccludedFormFactor(Pr, Nr, S);
  139.     else {
  140.         Frs = computeUnoccludedFormFactor(Pr, Nr, S);
  141.         if(Frs <= eps_F || quadArea(S) <= eps_A)
  142.             return(Frs * vis);
  143.         else {
  144.             splitQuad(S, dS);
  145.             Frs = 0.0;
  146.             for(j = 0; j < 4; j++)
  147.                 Frs += computeFormFactor(Pr, Nr, dS[j], Ns);
  148.         }
  149.     }
  150.  
  151.     return(Frs);
  152. }
  153.  
  154. /*
  155.     computeUnoccludedFormFactor--takes as input the descriptions of
  156.     a source patch and a receiving point and returns the unoccluded
  157.     form-factor from a differential area centered around the
  158.     receiving point to the finite area source.  The source is
  159.     described as a quadrilateral 'S' and its UNIT normal 'Ns'.
  160.     The receiving point is described by its location 'Pr' and its
  161.     UNIT normal 'Nr'. The form-factor is computed analytically.
  162. */
  163. static float computeUnoccludedFormFactor(Pr, Nr, S)
  164. Point Pr;
  165. Normal Nr;
  166. Quad S;
  167. {
  168.     int i;
  169.     float f, c;
  170.     Vector s, t, sxt, s0;
  171.  
  172.     f = 0.0;
  173.     VSUB(s, S[3], Pr);
  174.     c = 1.0 / VNORM(s);
  175.     VSCALE(s, c);
  176.     VCOPY(s0, s);
  177.     for(i = 0; i < 4; i++) {
  178.         if(i < 3) {
  179.             VSUB(t, S[i], Pr);
  180.             c = 1.0 / VNORM(t);
  181.             VSCALE(t, c);
  182.         } else
  183.             VCOPY(t, s0);
  184.         VCROSS(sxt, s, t);
  185.         c = 1.0 / VNORM(sxt);
  186.         VSCALE(sxt, c);
  187.         f -= acos(VDOT(s, t)) * VDOT(sxt, Nr);
  188.         VCOPY(s, t);
  189.     }
  190.     
  191.     return(f / (2.0 * M_PI));
  192. }
  193.  
  194. /*
  195.     splitQuad--takes as input a quadrilateral 'Q', splits it into
  196.     four equal quadrilaterals, and returns the result in array 'dQ'.
  197. */
  198. static void splitQuad(Q, dQ)
  199. Quad Q, dQ[4];
  200. {
  201.     int i, j;
  202.     Point center, midpt[4];
  203.  
  204.     /* compute the center of 'Q' and the midpoint of its edges */
  205.     VZERO(center);
  206.     for(i = 0; i < 4; i++) {
  207.         j = (i + 1) % 4;
  208.         VINCR(center, Q[i]);
  209.         VAVG2(midpt[i], Q[i], Q[j]);
  210.     }
  211.     VSCALE(center, 0.25);
  212.  
  213.     /* initialize the four new quadrilaterals */
  214.     VCOPY(dQ[0][0], Q[0]);
  215.     VCOPY(dQ[0][1], midpt[0]);
  216.     VCOPY(dQ[0][2], center);
  217.     VCOPY(dQ[0][3], midpt[3]);
  218.  
  219.     VCOPY(dQ[1][1], Q[1]);
  220.     VCOPY(dQ[1][2], midpt[1]);
  221.     VCOPY(dQ[1][3], center);
  222.     VCOPY(dQ[1][0], midpt[0]);
  223.  
  224.     VCOPY(dQ[2][2], Q[2]);
  225.     VCOPY(dQ[2][3], midpt[2]);
  226.     VCOPY(dQ[2][0], center);
  227.     VCOPY(dQ[2][1], midpt[1]);
  228.  
  229.     VCOPY(dQ[3][3], Q[3]);
  230.     VCOPY(dQ[3][0], midpt[3]);
  231.     VCOPY(dQ[3][1], center);
  232.     VCOPY(dQ[3][2], midpt[2]);
  233. }
  234.  
  235. /*
  236.     quadArea--takes as input a quadrilateral 'Q' and returns its
  237.     area.  The area of the quadrilateral is computed as the sum of
  238.     the areas of the two triangles obtained by splitting the
  239.     quadrilateral along one of its diagonals.
  240. */
  241. static float quadArea(Q)
  242. Quad Q;
  243. {
  244.     float area;
  245.     Vector u, v, uxv;
  246.  
  247.     VSUB(u, Q[2], Q[1]);
  248.     VSUB(v, Q[0], Q[1]);
  249.     VCROSS(uxv, u, v);
  250.     area = VNORM(uxv);
  251.  
  252.     VSUB(u, Q[0], Q[3]);
  253.     VSUB(v, Q[2], Q[3]);
  254.     VCROSS(uxv, u, v);
  255.     area += VNORM(uxv);
  256.  
  257.     return area * 0.5;
  258. }
  259.  
  260.